perm filename FLOAT.LSP[TIM,LSP] blob
sn#720407 filedate 1983-07-18 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00008 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00002 00002 Load MACHAR routine
C00011 00003 Random number generators
C00012 00004 Square root program
C00016 00005 Square-root Testing
C00022 00006 ARCTAN/ARCTAN2
C00028 00007 ARCTAN Test Routine
C00039 00008 Testing
C00040 ENDMK
C⊗;
;;; Load MACHAR routine
(eval-when (compile load) (fasload machar))
(declare (special *ibeta* *it* *irnd* *ngrd* *machep* *epsneg* *maxexp*
*negep* *iexp* *minexp* *eps* *xmin* *xmax*))
(declare (flonum *xmax* *xmin* *epsneg*)
(fixnum *ibeta* *it* *irnd* *machep* *minexp* *iexp*
*negep* *ngrd* *maxexp*))
(declare (*expr machar))
;;; Random number generators
(declare (special *iy* *j*) (fixnum *iy* *j*)
(flonum (ran fixnum))
(flonum (randl flonum)))
(setq *iy* 100001.)
(defun ran (k)
(setq *j* k)
(setq *iy* (* *iy* 125.))
(setq *iy* (- *iy* (* (// *iy* 2796203.) 2796203.)))
(//$ (float *iy*) 2796203.0))
;;; Get value of E in here
(defun randl (x)
(exp (*$ x (ran 1))))
;;; Square root program
;;; intxp(x) returns the integer representation of the
;;; exponent in the normalized representation of the floating-point
;;; number x. intxp(3.0)= 2
;;; adx(x,n) augments the integer exponent in the floating-point
;;; representation of X by N, thus scaling X by the Nth power of the
;;; radix. adx(1.0,2) = 4.0
;;; setxp(x,n) returns the floating-point representation of
;;; a number whose significand is the significand of X, and whose exponent
;;; is N. setxp(1.0,3)=4
(declare (fixnum (intxp flonum))
(flonum (adx flonum fixnum))
(flonum (setxp flonum fixnum))
(*expr adx setxp intxp))
(defun square-root (x)
(declare (flonum x f y)
(fixnum n i))
(cond ((zerop x) 0.0)
((< x 0.0)
(terpri)
(princ "Square-root of a negative number")
(terpri)
0.0)
(t (let ((n (intxp x))
(f (setxp x 0)))
(let ((y (+$ .41731 (*$ .59016 f))))
(do ((i 3 (1- i))
(y y (*$ .5 (+$ y (//$ f y)))))
((zerop i)
(cond ((oddp n)
(setq n (1+ n))
(setq y (*$ y
.70710678118654752440))))
(adx y (// n 2)))))))))
;;; Lap examples of how to implement the functions above:
(lap intxp subr)
(args intxp (nil . 1))
(push p (% 0 0 fix1)) ;skip this entry if declared
(move t 0 a) ;gets the machine representation of x
(ldb tt bpt) ;extracts the exponent (bits 1-8 on a PDP-10)
(subi tt #o200) ;the exponent is stored excess #o200
(popj p) ;fxcons it, maybe
(entry adx subr)
(args adx (nil . 2))
(push p (% 0 0 float1))
(move tt 0 a) ;get x
(ldb t bptt) ;get exponent
(add t 0 b) ;add n
(dpb t bptt) ;put the new exponent in place
(popj p) ;float CONS it (maybe)
(entry setxp subr)
(args setxp (nil . 2))
(push p (% 0 0 float1)) ;skip this if doing declared call
(move tt 0 a) ;get x
(move t 0 b) ;get n
(addi t #o200) ;make up excess #o200 exponent
(dpb t bptt) ;install the exponent
(popj p) ;float CONS (maybe)
bpt (331000←22 0 t);byte pointers
bptt (331000←22 0 tt)
()
;;; Square-root Testing
(declare (special *results*))
(setq *results* ())
(defun sqrt-test ()
(declare (flonum beta sqbeta albeta ait a b xn c x1 r6 r7 x y z w)
(fixnum n i k1 k2 k3))
(machar)
(let* ((beta (float *ibeta*))
(sqbeta (square-root beta))
(albeta (log beta))(ait (float *it*))
(a (//$ 1.0 sqbeta))(b 1.0)
(n 8000.)(xn (float n))
(c 0.0)(k1 0)(k3 0)(x1 0.0)(r6 0.0)(r7 0.0)
(x 0.0)(y 0.0)(z 0.0)(w 0.0)(k2 0))
(do ((j 1 (1+ j)))
((= j 3))
(setq c (log (//$ b a)))
(setq k1 0 k3 0 x1 0.0 r6 0.0 r7 0.0)
(do ((i 1 (1+ i)))
((> i n))
(setq x (*$ a (randl c)))
(setq y (*$ x x))
(setq z (square-root y))
(setq w (//$ (-$ z x) x))
(cond ((> w 0.0) (setq k1 (1+ k1))))
(cond ((< w 0.0) (setq k3 (1+ k3))))
(setq w (abs w))
(cond ((< r6 w)
(setq r6 w)
(setq x1 x)))
(setq r7 (+$ r7 (*$ w w))))
(setq k2 (- (- n k1) k3))
(setq r7 (square-root (//$ r7 xn)))
(push '(test of sqrt(x * x) - x) *results*)
(push
`(,n random arguments were tested in the interval
(,a ,b)) *results*)
(push
`(sqrt(x) was larger ,k1 times) *results*)
(push `(it agreed ,k2 times) *results*)
(push `(it was smaller ,k3 times) *results*)
(push `(there are ,*it* base ,*ibeta* significant digits in
a floating-point number) *results*)
(setq w -999.0)
(cond ((not (zerop r6))
(setq w (//$ (log (abs r6)) albeta))))
(push `(the maximum relative error of ,r6
= ,*ibeta* ↑ ,w occurred for x =
,x1) *results*)
(setq w (max (+$ ait w) 0.0))
(push `(the estimated loss of base ,*ibeta*
significant digits is
,w) *results*)
(setq w -999.0)
(cond ((not (zerop r7))
(setq w (//$ (log (abs r7)) albeta))))
(push `(the root mean square relative
error was ,r7 = ,*ibeta* ↑ ,w) *results*)
(setq w (max (+$ ait w) 0.0))
(push `(the estimated loss of base ,*ibeta*
significant digits is
,w) *results*)
(setq a 1.0)
(setq b sqbeta))
(push `(test of special arguments) *results*)
(setq x *xmin*)
(setq y (square-root x))
(push
`(sqrt(*xmin*) = sqrt(,x) = ,y) *results*)
(setq x (-$ 1.0 *epsneg*))
(setq y (square-root x))
(push `(sqrt (1.0 - *epsneg*) = sqrt(1.0 - ,*epsneg*) = ,y) *results*)
(setq x 1.0)
(setq y (square-root x))
(push `(sqrt(1.0) = ,y) *results*)
(setq x (+$ 1.0 *eps*))
(setq y (square-root x))
(push `(sqrt (1.0 + *eps*) = sqrt (1.0 + ,*eps*) = ,y) *results*)
(setq x *xmax*)
(setq y (square-root x))
(push `(sqrt(*xmax*) = sqrt(,*xmax*) = ,y) *results*)
(push '(test of error returns) *results*)
(setq x 0.0)
(push `(sqrt will be called with an argument of ,x this
should not trigger an error) *results*)
(setq y (square-root x))
(push `(sqrt returned the value ,y) *results*)
(setq x -1.0)
(push `(sqrt will be called with an argument of ,x this
should trigger an error) *results*)
(setq y (square-root x))
(push `(sqrt returned the value ,y) *results*)
(push '(this concludes the tests) *results*)
t))
;;; ARCTAN/ARCTAN2
;;; Set EPS to *ibeta*↑(*it*/2)
(eval-when (eval compile)
(machar)
(setq eps (expt (float *ibeta*) (//$ (float *it*) -2.0)))
;;;Define R depending on the number of bits in the significand.
(cond ((or (and (= 2 *ibeta*)
(< *it* 25.))
(and (= 10. *ibeta*)
(< *it* 8.)))
(defmacro R (x)
`(let ((q (*$ (+$ (*$ -0.5090958253e-1 ,x)
-0.4708325141) ,x))
(r (+$ ,x
0.1412500740e1)))
(//$ q r))))
((or (and (= 2 *ibeta*)
(< *it* 33.))
(and (= 10. *ibeta*)
(< *it* 11.)))
(defmacro R (x)
`(let ((q
(*$ (+$ (*$ -0.720026848898 ,x)
-0.144008344874e1)
,x))
(r (+$ (*$ (+$ ,x
0.475222584599e1)
,x)
0.432025038919e1)))
(//$ q r))))
((or (and (= 2 *ibeta*)
(< *it* 51.))
(and (= 10. *ibeta*)
(< *it* 16.)))
(defmacro R (x)
`(let ((q
(*$ (+$ (*$ (+$
(*$ -0.794391295408336251 ,x)
-0.427444985367930329e1)
,x)
-0.427432672026241096e1)
,x))
(r (+$ (*$ (+$ (*$
(+$ ,x 0.919789364835039806e1)
,x)
0.205171376564218456e2)
,x)
0.128229801607919841e2)))
(//$ q r))))
((or (and (= 2 *ibeta*)
(< *it* 61.))
(and (= 10. *ibeta*)
(< *it* 19.)))
(defmacro R (x)
`(let ((q
(*$ (+$ (*$ (+$ (*$ (+$ (*$
-0.83758299368150059274 ,x)
-0.84946240351320683534e1)
,x)
-0.20505855195861651981e2)
,x)
-0.13688768894191926929e2)
,x))
(r (+$ (*$ (+$ (*$ (+$ ,x 0.15024001160028576121e2)
,x)
0.59578436142597344465e2)
,x)
0.86157349597130242515e2)))
(//$ q r))))))
(defmacro arcfinishb ()
'(cond ((< u 0.0)
(setq result
(-$ 3.14159265358979323846 result))
(arcfinishc))
(t (arcfinishc))))
(defmacro arcfinishc ()
'(cond ((< v 0.0)
(minus result))
(t result)))
(defun arctan (x)
(declare (flonum x))
(arctanf x () 0.0 0.0))
(defun arctan2 (v u)
(declare (flonum u v result)
(fixnum expdiff))
(cond ((zerop u)
(cond ((zerop v)
(terpri)
(princ "arctan2 called with u = 0.0 and v = 0.0")
(terpri) 0.0)
(t (let ((result
1.57079632679489661923))
(arcfinishc)))))
(t (let ((expdiff (- (intxp u)(intxp v))))
(cond
((not (< expdiff #.(- *maxexp* 3))) ;overflow?
(let ((result
1.57079632679489661923))
(arcfinishc)))
((not (< #.(+ *minexp* 3) expdiff)) ;underflow?
(let ((result 0.0))
(arcfinishb)))
(t (arctanf (//$ v u) t u v)))))))
(defun arctanf (x arctan2p U V)
(declare (flonum f g result u v x)
(fixnum n))
(let ((n 0) (result 0.0)(f (abs x)))
(cond ((> f 1.0)
(setq f (//$ 1.0 f))
(setq n 2)))
(cond ((> f 0.26794919243112270647)
(setq f (//$ (+$ (-$ (-$ (*$ 0.73205080756887729353 f)
0.5)
0.5) f)
(+$ 1.73205080756887729353 f)))
(setq n (1+ n))))
(cond ((< (abs f) #.eps)
(setq result f))
(t (let ((g (*$ f f)))
(setq result (+$ f (*$ f (r g)))))))
(cond ((> n 1) (setq result (minus result))))
(setq result (+$
(caseq n
(0 0.0)
(1 0.52359877559829887308)
(2 1.57079632679489661923)
(t 1.04719755119659774615))
result))
(cond (arctan2p
(arcfinishb))
(t (cond ((< x 0.0)
(minus result))
(t result))))))
;;; ARCTAN Test Routine
(defun arctan-test ()
(declare (flonum beta albeta ait a b betap del em expon ob32 ran r6 r7
xn w x xl xsq x1 y z zz)
(fixnum n i k1 k2 k3 ii i1 j))
(machar)
(let* ((beta (float *ibeta*))
(albeta (log beta))(ait (float *it*))
(a -0.0625)(b (minus a))(ob32 (*$ b 0.5))
(n 8000.)(xn (float n))(x1 0.0)(xl 0.0)(expon 0.0)
(r7 0.0)(r6 0.0)(x 0.0)(z 0.0)(xsq 0.0)(em 0.0)
(k1 0)(k3 0)(del 0.0)(sum 0.0)(zz 0.0)
(i1 0)(y 0.0)(w 0.0)(k2 0)(betap 0.0))
(do ((j 1 (1+ j)))
((> j 4))
(setq k1 0)(setq x1 0.0)(setq r7 0.0)
(setq k3 0)(setq r6 0.0)(setq del (//$ (-$ b a) xn))
(setq xl a)
(do ((i 1 (1+ i)))
((> i n))
(setq x (+$ (*$ del (ran i1))
xl))
(cond ((= j 2)
(setq x (*$ (-$ (+$ 1.0 (*$ x a)) 1.0) 16.0))))
(setq z (arctan x))
(cond ((= j 1)
(setq xsq (*$ x x))
(setq em 17.0)
(setq sum (//$ xsq em))
(do ((ii 1 (1+ ii)))
((> ii 7))
(setq em (-$ em 2.0))
(setq sum (*$ (-$ (//$ 1.0 em) sum) xsq)))
(setq sum (*$ (minus x) sum))
(setq zz (+$ x sum))
(cond ((zerop *irnd*) (setq zz (+$ zz (+$ sum sum))))))
(t (cond ((= j 2)
(setq y (-$ x .0625))
(setq y (//$ y (+$ 1.0 (*$ x a))))
(setq zz (+$ (-$ (arctan y)
8.11900402651526021e-5)
ob32))
(setq zz (+$ zz ob32)))
(t (setq z (+$ z z))
(setq
y (//$ x (*$ (+$ 0.5 (*$ x 0.5))
(+$ (-$ 0.5 x) 0.5))))
(setq zz (arctan y))))))
(setq w 1.0)
(cond ((not (zerop z))
(setq w (//$ (-$ z zz) z))))
(cond ((> w 0.0)
(setq k1 (1+ k1))))
(cond ((< w 0.0)
(setq k3 (1+ k3))))
(setq w (abs w))
(cond ((< r6 w)
(setq r6 w)
(setq x1 x)))
(setq r7 (+$ r7
(*$ w w)))
(setq xl (+$ xl del)))
(setq k2 (- (- n k3) k1))
(setq r7 (sqrt (//$ r7 xn)))
(cond ((= j 1)
(push `(Test of ARCTAN ( x ) vs truncated
taylor series)
*results*)))
(cond ((= j 2)
(push `(Test of ARCTAN ( x )
vs arctan ( 1 // 16.) +
arctan((x - 1 // 16.)//(1 + x // 16.)))
*results*)))
(cond ((= j 3)
(push `(test of 2 * arctan ( x ) vs arctan
(2x // (1 - x * x)))
*results*)))
(push `(,n random arguments were tested from the
interval ( ,a ,b))
*results*)
(push `(arctan ( x ) was larger ,k1 times)
*results*)
(push `(it agreed ,k2 times) *results*)
(push `(it was smaller ,k3 times) *results*)
(push `(there are ,*it* significant base ,*ibeta* digits in
a floating-point number) *results*)
(setq w -999.0)
(cond ((not (zerop r6))
(setq w (//$ (log (abs r6)) albeta))))
(push `(the maximum relative error of ,r6 = ,*ibeta*
↑ ,w occurred for x = ,x1)
*results*)
(setq w (max (+$ ait w) 0.0))
(push `(the estimated loss of base ,*ibeta* significant digits
is ,w) *results*)
(setq w -999.0)
(cond ((not (zerop r7))
(setq w (//$ (log (abs r7)) albeta))))
(push `(the root mean square relative error was ,r7
= ,*ibeta* ↑ ,w) *results*)
(setq w (max (+$ ait w) 0.0))
(push `(the estimated loss of base ,*ibeta* significant digits
is ,w) *results*)
(setq a b)
(cond ((= j 1) (setq b (-$ 2.0 (sqrt 3.0)))))
(cond ((= j 2) (setq b (-$ (sqrt 2.0) 1.0))))
(cond ((= j 3) (setq b 1.0))))
(push `(special tests) *results*)
(push `(the identity: arctan (-x) = -arctan (x) will be tested)
*results*)
(push `(x : f (x) + f (-x)) *results*)
(setq a 5.0)
(do ((i 1 (1+ i)))
((> i 5))
(setq x (*$ (ran i1) a))
(setq z (+$ (arctan x)
(arctan (minus x))))
(push `(,x : ,z) *results*))
(push `(the identity arctan (x) = x for x small will be tested)
*results*)
(push `(x : x - f ( x)) *results*)
(setq betap (expt beta *it*))
(setq x (//$ (ran i1) betap))
(do ((i 1 (1+ i)))
((> i 5))
(setq z (-$ x (arctan x)))
(push `(,x : ,z) *results*)
(setq x (//$ x beta)))
(push `(the identity arctan (x // y) = arctan2 (x y) will
be tested) *results*)
(push `(the first column of results should be 0 and the
second should be +-π) *results*)
(push `(x : y : f1 ( x // y) - f2 (x y) : f1 (x // y) - f2 (x // -y))
*results*)
(setq a -2.0)
(setq b 4.0)
(do ((i 1 (1+ i)))
((> i 5))
(setq x (+$ (*$ (ran i1) b) a))
(setq y (ran i1))
(setq w (minus y))
(setq z (-$ (arctan (//$ x y)) (arctan2 x y)))
(setq zz (-$ (arctan (//$ x w)) (arctan2 x w)))
(push `(,x : ,y : ,z : ,zz) *results*))
(push `(test of very small argument) *results*)
(setq expon (*$ (float *minexp*) 0.75))
(setq x (expt beta expon))
(setq y (arctan x))
(push `(arctan ( ,x ) = ,y) *results*)
(push `(test of error returns) *results*)
(push `(arctan will be called with the argument ,*xmax*) *results*)
(push `(this should not trigger an error message) *results*)
(setq z (arctan *xmax*))
(push `(arctan ( ,*xmax* ) = ,z) *results*)
(setq x 1.0 y 0.0)
(push `(arctan2 will be called with the arguments ,x ,y) *results*)
(push `(this should not trigger an error message) *results*)
(setq z (arctan2 x y))
(push `(arctan2 ( ,x ,y) = ,z) *results*)
(push `(arctan2 will be called with the arguments ,*xmin*
,*xmax*) *results*)
(push `(this should not trigger an error message) *results*)
(setq z (arctan2 *xmin* *xmax*))
(push `(arctan2 ( ,*xmin* ,*xmax*) = ,z) *results*)
(push `(arctan2 will be called with the arguments ,*xmax*
,*xmin*) *results*)
(push `(this should not trigger an error message) *results*)
(setq z (arctan2 *xmax* *xmin*))
(push `(arctan2 ( ,*xmax* ,*xmin*) = ,z) *results*)
(setq x 0.0)
(push `(arctan2 will be called with the arguments ,x ,y)
*results*)
(push `(this should trigger an error message) *results*)
(setq z (arctan2 x y))
(push `(arctan2 ( ,x ,y) = ,z) *results*)
(push `(This concludes the tests) *results*)
t))
;;; Testing
(defun show-results ()
(mapc #'print (reverse *results*)) t)
(include "timer.lsp[tim,lsp]")
(timer timit1 (progn (sqrt-test)(arctan-test)))
(defun timit ()
(timit1)
(show-results))